A microscopic 2D lattice model of dimer granular compaction with friction 
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We study by Monte Carlo simulation the compaction dynamics of hard dimers in 2D under 
the action of gravity, subjected to vertical and horizontal shaking, considering also the case in 
which a friction force acts for horizontal displacements of the dimers. These forces are modeled by 
introducing effective probabilities for all kinds of moves of the particles. We analyze the dynamics 
for different values of the time r during which the shaking is applied to the system and for different 
intensities of the forces. It turns out that the density evolution in time follows a stretched exponential 
behavior if r is not very large, while a power law tail develops for larger values of r. Moreover, in 
the absence of friction, a critical value r* exists which signals the crossover between two different 
regimes: for r < r* the asymptotic density scales with a power law of r, while for t > r* it reaches 
logarithmically a maximal saturation value. Such behavior smears out when a finite friction force 
is present. In this situation the dynamics is slower and lower asymptotic densities are attained. 
In particular, for significant friction forces, the final density decreases linearly with the friction 
coefficient. We also compare the frictionless single tap dynamics to the sequential tapping dynamics, 
observing in the latter case an inverse logarithmic behavior of the density evolution, as found in the 
experiments. 
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I. INTRODUCTION 

The process of compaction in granular media attracts 
a great deal of attention in the scientific community, 
in particular amongst physicists and chemical engineers. 
In fact it displays features that are general enough to 
make it suitable for investigation through models that 
are relatively simple, at least for a non equilibrium pro- 
cess. Granular compaction therefore appears appealing 
for testing new promising and unifying ideas in the field 
of disordered systems flfl]. 

Grains poured into a vessel fill it with loose arrange- 
ments and relatively low densities. The action of external 
perturbations, like shaking and tapping, in the presence 
of an external driving, like the gravity field, leads to a 
very slow increase of density through a rich phenomenol- 
ogy, displaying different regimes and both reversible and 
irreversible dynamical phases |5|-|l8f|. A number of dif- 
ferent lattice models have been recently investigated in 
order to unravel the microscopic mechanisms producing 
such behaviors, and to clarify to which extent they may 
be considered general, also in relation to other disordered 



systems like glasses |l9|-|27|]. 

In this work we introduce a lattice gas model for dimer 
compaction which includes, besides the vertical, also hor- 
izontal shaking and friction. Friction, although playing a 
fundamental role in real granular systems, has not been 
introduced so far in models for compaction. Horizon- 
tal shaking has been considered only experimentally for 
non Brownian spheres where it seems to favor local dense 
crystalline order 28 29[|. Highly anisotropic granulars, 
namely rods, have been recently investigated under ver- 
tical shaking. In contrast to sphere packings which tend 
to end up in disordered configurations, evidence is found 
that the particle anisotropy drives ordering |50| . We try 
to incorporate in our model, albeit in a simplified man- 
ner, all these aspects, namely particle anisotropy, hori- 
zontal shaking and frictional dynamics, which we find to 
affect significantly the compaction process. 

In particular, our model is based on the diffusional 
dynamics of dimers, considered as rigid and non over- 
lapping particles which occupy two lattice sites on a 
square lattice. The model has no quenched disorder, 
but for this type of dynamics glass-like properties may 
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be expected [pi[ , because of the onset of geometrically 
frustrated configurations. Since the number of states 
with the closest packed density is exponentially large 
in the lattice size ^||,^3), this model is intrinsically dif- 
ferent from the other lattice models recently investi- 
gated [pj|-p5| and is suitable to understand how horizon- 
tal shaking and friction may influence the general char- 
acteristics attributed to the compaction process. 

In the next section we describe in detail the model and 
the Monte Carlo method that we used. In Sec. Ill 



discuss the results for the compaction dynamics without 
friction and illustrate the evolution of the packing struc- 
In Sec. IV we introduce the effect of friction and 

The final 



ture. 



show its crucial role in the dynamical behavior, 
section is devoted to the conclusions. 



II. MODEL 



We consider a square lattice of N = L x L sites, with 
lateral periodic boundary conditions, an open boundary 
at the top and a rigid wall at the bottom. Elongated par- 
ticles occupying two consecutive lattice sites (dimers) are 
inserted from the top one at the time with random (hor- 
izontal or vertical) orientation, letting them fall down 
keeping their orientation fixed, until reaching a stable 
position. 



FIG. 1. Schematic picture of the lattice model studied. 
The arrangement of the particles is the result of the initial 
sample preparation. Note that lateral periodic boundary con- 
ditions have been imposed. 



In this way we are able to prepare the system in an ini- 
tial state which is saturated, i.e. no more particles can be 
put in, and with statistically reproducible density po, i.e. 
characterized by a precise mean value, po ~ 0.587, cor- 
responding to a random loose packing. In the present 
model, particles are subjected only to geometrical in- 
teractions, and the non-overlapping condition produces 
strong constraints on their relative positions (see Fig. [l]). 
The aim of this paper is to study the system in presence 
of gravity and external vibrations, with the possibility of 
taking into account also a friction-like force for horizontal 
displacements. For this purpose we perform Monte Carlo 
(MC) simulations introducing a random diffusive dynam- 
ics in our model, which mimics the aforementioned forces 
and preserves the geometrical constraints. At this stage, 
we consider separately horizontal and vertical dimers and 
keep their orientation fixed without allowing rotations. 
We plan to introduce this feature, which would make the 
model more realistic, in the future. We consider in de- 
tail a single tap applied to the system for a fixed time t, 
comparing it with the multiple tapping for selected cases 
(see |l9j). We first describe the dynamics without fric- 
tion. The dynamics can be divided into two stages: (i) for 
t < t particles can move horizontally (left or right with 
probability p/i/2) and vertically (upwards with probabil- 
ity p up or downwards with probability Pdown)', (ii) for 
t > t only downward movement of particles is possible 
(i.e. p up = p h = and pdown = 1)- We note that t _1 
can be thought of as the quench rate of an initially an- 
nealed system. The probabilities for the different moves 
are normalized: p h + p up + Pdown = 1- Physically Pdown 
corresponds to the action of gravity, while ph and p up 
represent respectively the horizontal and vertical compo- 
nent of the external tap. The functional form chosen for 
the time dependence of ph and p up is pj| 



Ph(t)=p h (l-t/r)9(r-t), 



(1) 



where 0(x) is the Heaviside function (the same expres- 
sion holds for p U p(t), with p up instead of pi). Pdown(t) is 
determined from the normalization condition: 



Pdown{t) = 1 -Pup(t) ~Ph{t) 



- r>° 
"down 



(p°u P +PhWr 
(2) 

where p a down = 1 - (p° up + Since p do wn{t) > Vf 

it must be p up + Pdown — 1- If not stated explicitly 
otherwise, we assume Pd own — Pup> so Ph results to be 
pi = 1 — 2p U pi that is only one independent input pa- 
rameter is needed. In each MC move one extracts a par- 
ticle with uniform probability, chooses a move for this 
particle according to the values of p%, p up and pdown, 
and performs the move if all the geometrical constraints 
are satisfied. One MC time step (MCS) corresponds to 
N attempted moves. In the following, time is always 
given in units of MCS. To save CPU time we used an 
algorithm in which the attempted moves are always suc- 
cessful, and consistently updated time through proba- 
bilistic arguments (for the details of our computation see 
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Ref. |34|] on a reaction-diffusion model for dimers). We 
performed our MC simulations on a lattice with L = 100, 
for which we checked that finite size effects are negligible. 

As a next step we introduce friction between adjacent 
horizontal layers in the material in the following simpli- 
fied way: the effective probability of making a horizontal 
move is set to 

Ph f (*) = Ph(t) ~ ^(Pdown(t) -Pup(t)) = 

= \p° h - M (l -p° h - 2p° up ) - \p° h + ix(p° h + 2p° up )}t/T]0(t o ( t x) 

(3) 



where 



*oO) 



Pi 



fx(l 



Pi 



p° h + fi(p° h + 2 Puf 



(4) 



p acts as a friction coefficient (p > 0) and the "friction 
force" is proportional to the load, represented by the net 
vertical force, reducing the probability of the move. This 
assumption implies that the frictional force occurs only 
for horizontal moves, independently of whether an un- 
derlying dimer (i.e. a dimer in the adjacent lower row) is 
present. We are aware that this is a simplified approxi- 
mation, but it is meant as a first step in describing this 
complex process. To ensure p'ff* > at t = 0, p has to 
satisfy the condition 



< 



Pi 



2p2p)' 



In addition, for the friction force to give a negative con- 
tribution to ph, one must have pdown(t) > p U p(t) Vf (this 
is automatically satisfied with our choice p down = P up )- 
Note that the normalization of the frictionless probabili- 
ties is still p up + pdown +Ph = 1 ■ I 11 other terms we have 
introduced a "sticking" probability (probability that the 
particle does not move) p s , given by p s = p{pdown ~Pup)- 
If we put p = we recover the frictionless case. 



III. COMPACTION DYNAMICS WITHOUT 
FRICTION 

We have analyzed the time behavior of the mean den- 
sity of the system p(t) (number of occupied sites normal- 
ized to the total number of sites) by measuring it in the 
lower 30% of the box at time t. In this way we are sure 
to measure density at the bulk, since fluctuations can in- 
crease significantly in the proximity of the open bound- 
ary, resulting in a quite complicated compaction behavior 
in this region (2^| . Since the statistical fluctuations of the 
density can be quite relevant we have performed several 
MC realizations (up to 500) of the process in order to 
obtain reliable results. After preparing the system in its 
low density configuration (p = p(t = 0) ~ 0.587), the 



diffusive dynamics starts as described in the previous sec- 
tion, p^ is chosen as input parameter and consequently 
Pup = Pdown = (1 — Ph)/2- We stop the simulations when 
a steady asymptotic value p^, of the density is reached. 
We have studied the density evolution for different values 
of jP h and of the shaking time r, in order to find out if 
some kind of scaling law, describing the behavior of the 
final density as a function of these parameters, exists. 
As far as we know, no systematic study concerning this 
point has been performed. 

- ^)The time behavior of the density for several choices of 
r is illustrated in Fig. ||, where a relatively low value of 
p\ has been used (p h = 0.1). 
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FIG. 2. Temporal behavior of density in the frictionless 
case 

for different values of the shaking time r (from left to right 

t = 10, 30, 50, 80, 103, 125, 150, 210, 280, 400, 800, 2000, 4000, 10000). 

The horizontal shaking amplitude is p° = 0.1. 

We have observed that the dynamics gets drastically 
slower when r increases; actually, because of our shaking 
rule, Eq. (|l|), when a long tap is applied the density does 
not significantly change in the first steps of the process, 
but the evolution takes place on a much wider time scale 
and finally a larger value of the asymptotic density is 
achieved. We have tried to fit different functional forms 
to the data in Fig. |^, looking in particular at the re- 
laxation functions proposed in the seminal experimental 
paper of Ref. ||. It is clear that the observed dynami- 
cal behavior is very complex and is not compatible with a 
single relaxation time, i.e. a simple exponential. We have 
found that the commonly claimed inverse logarithmic re- 
laxation does not hold for our system in the single tap 
case. Instead, the most suitable functional form for our 
data is a stretched exponential, as proposed by Nicodemi 
et al. |l9[ ] for the single tapping dynamics: 

p{t- T,p° h ) = Poc(r,p° h ) - C{ Poo {t,pD - po) exp[-((f + t Q )/r f] 

(5) 
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where the fit parameters C, to, t$ and (3 depend in prin- 
ciple both on r and Actually, it turns from our fits 
that C ~ 1.1 and /3 ~ 10 are almost independent of r 
and In Fig. ||(a) we show the density evolution for a 
small value of r. 




10' 10' t(MCS) 10° 10- 

FIG. 3. Density relaxation in the frictionless case for 
t = 80 (a) and r = 4000 (b). The solid lines are the re- 
sult of the MC simulation, the dashed lines are stretched ex- 
ponential fits using Eq. (|5j) for p(t) > 0.6, with parameters 
C ~ 1.1, t ^ 170, r ~ 250, (3 ~ 9.8 for r = 80 (a), and 
C ~ 1.1, to ^ 2900, r ~ 3600, f3 ~ 9.7 for r = 4000 (b). 
The dot-dashed line in (b) is a fit for t > 1000 according to 
Eq. (|) with parameters B ~ 0.018, r ~ 330 and a ~ 4.95. 
The horizontal shaking amplitude is p° = 0.1. 

As it can be seen, the intermediate-long time behav- 
ior can be accurately fitted by Eq. (Q) (all fits have been 
performed for p > 0.6). As t increases, however, Eq. (||) 
fails to reproduce the long time regime. In fact, a power 
law tail develops, indicating that the compaction process 
slows down at long times for large shaking duration. The 
fitting function that we have used for the long time be- 
havior is 



o(t-Tv°)-o (tv°) P^ p h)-P* 



(6) 



This is illustrated in Fig. |(b) for r = 4000. Fit (§) has 
been applied for r > 1000. The exponent a of this power 
law can be considered independent of r (a ~ 5). How- 
ever the stretched exponential function still describes the 
density behavior in the intermediate time regime well. 

We have also studied the density relaxation for differ- 
ent values of p° at fixed r (r = 103). The corresponding 
plots are shown in Fig. [|, where it is clear that the role 
of the horizontal shaking amplitude is crucial in deter- 
mining the asymptotic density. 
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FIG. 4. Temporal behavior of density in the frictionless 
case for different values of the horizontal shaking amplitude 
pi (from bottom to top p° h = 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1). The 
shaking time is r = 103. 



Moreover we have considered the effect of different ra- 
tios Pup/Pdown on ^ ne density evolution, which we show 
in Fig. H When p^ p > p down we observe a decompaction 
at the beginning of the process, which becomes more pro- 
nounced if the vertical tap is stronger. 



o, o , 

Pup 'Pdown ~ 1 

Pup 'Pdown 
0,„ 
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FIG. 5. Temporal behavior of density in the frictionless 
case for different vafues of the ratio p° p /p° oron : pZ P /Pdown = 1 
(solid line), pl p /p° down = 2 (dashed line) and pl p /p° down = 3 
(dotted line). The shaking time is r = 80 and the horizontal 
shaking amplitude is p° = 0.1. 



This reflects the fact that, in the initial transient, den- 
sity decreases at the bottom of the container for a strong 
tap. Such a decompaction upon vertical acceleration in- 
crease seems to be a genuine feature of two-dimensional 
systems, as claimed in | p5[ . The saturation density is also 
sensitive to variations of p^ p and decreases for increas- 
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ing Pup/Pdown- However, the scaling behavior does not t=o(MCS) .. :=3775 (Mcs; 

change and the density evolution can still be described 
by Eq. ||. but with a smaller value of 3 with respect to 
tlie case p° up = p down (0 ~ 9.8 for p° up /p down = 1,0* 7.6 
for vl P IP down = 2 and * 6.5 for p° up /p down = 3). 

In order to gain a better understanding of the degree 
of compaction of the system, we investigated the depen- 
dence of the saturation density /< x on r and p\. The 

data for p^ vs. t are plotted in Fig. ^(a). ______________ 

1=1 026.7 (MCS) t=4032.7 (MCS) [end] 
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FIG. 7. Snapshots of one particular realization of the pack- 
ing at different times for r = 4000. The horizontal shaking 
amplitude is p° = 0.1. 



FIG. 6. Saturation density p x in the frictionless case as 
a function of the shaking time r (a) and the ratio pl/Pdown 
(b). The squares are the numerical data, while the solid line 
and dotted line in (a) are respectively the fits for r < r* and 
r > r* given in Eq. |7J The horizontal shaking amplitude in 
(a) is p° h = 0.1 and the shaking time in (b) is r = 103. 

A crossover from a power law scaling to a logarithmic 
behavior is observed at a characteristic value t* ~ 200: 

Poo(t) = {t/ti) s forr<T* 

Poo(t) = p m - B/\n(r/T r ) for r > r* (7) 

where n ~ 4.5 • 10 3 , 6 ~ 0.073, p m ~ 1, B ~ 0.6 and 
r r ~ 18. This means that for large values of r the pro- 
cess gets slower and slower, as observed before, eventually 
reaching a final value of the density p m ~ 1 for r — > oo 
(which corresponds to about the maximal density for the 
system p4j). The variation of p^ with respect to the 
ratio PhlPdown i s illustrated in Fig. ||(b). The increase of 
Poo is more pronounced for low values of p® , while it tends 
to saturate afterwards. This compares qualitatively well 
to the experimental results for a horizontally shaken box 
filled with beads pq] , where it was found that for low 
filling rates the packing crystallizes upon increase of the 
adimensional parameter T — Auj 2 /g (where A and u> are 
the amplitude and pulsation of the vibration and g is the 
gravitational acceleration) , which roughly corresponds to 

Phi Pdown- 

In Fig. M we show some snapshots of one particular 
realization of the compaction process for r — 4000. 



We see how the packing evolves from a highly disor- 
dered state at t = to a polycrystalline state, made of 
many crystalline ordered domains, at the end of the com- 
paction. As time increases the order proceeds from the 
bottom to the top of the packing while empty spaces re- 
duce in size. The final state is therefore characterized 
by the presence of large clusters of horizontal and ver- 
tical dimers. In order to demonstrate that the shaking 
amplitude p {) h and time r play a significant role, we have 
reported in Fig. || a comparison between the final configu- 
rations for p° h = 0.1 and p° h = 0.8 at fixed t — 103 (upper 
part), and for r = 80 and t = 4000 at fixed p° = 0.1 
(lower part). 
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FIG. 8. Snapshots of the final configuration of the packing 
for two different pi and fixed r (upper part) and two different 
t and fixed p\ (lower part) . 

We notice that an intense or/and long tap is effective 
in locally removing frozen configurations through collec- 
tive rearrangements of particles, thus letting the defects 
migrate towards the free surface of the packing. 

We have also investigated in detail the behavior of the 
parameters of the stretched exponential fit (Eq. (||)) to 
and tq, which is reported in Fig. ^. 



same value determining the change in the behavior of 
Poo vs. r described in Eq. (Q). A similar transition was 
also signaled by Nicodemi et al. jl9|. The parameters 
of the stretched exponential relaxation Eq. (||) change 
according to: 

t (r) = r/ti + t< for t < r* 

*o(t) = (r/hy + a for t>t* (8) 

to(t) = t/t\ + t< for T < T* 
to(t) = (t/t 2 ) 7 + b for t>t* (9) 

where h ~ 1.40, i< ~ 130, t 2 ~ 7 • 10" 5 , 7 ~ 0.47, 
a ~ -700, n ~ 0.76, t< ~ 150, r 2 ~ 10~ 4 , b ~ -700. 
It is remarkable that both to and To vary linearly with r 
for t < t* and both follow a power law with the same 
exponent 7 for r > r*. This could mean that they are 
both determined by the same relaxation mechanism as 
r changes (i.e. if we define to = to — t<, fo = To — t<, 
?q = to — a and fo = To — b, we see from Eqs. (0) and (g) 
that to/fo and to/ fa are constants independent of t). 

We note that the critical shaking time t* physically 
signals the crossover between two different regimes. In 
fact, a careful inspection of Fig. § shows that, for r < t* , 
compaction occurs after the shaking (when = p up = 
and Pdown — 1) and is mainly driven by gravity. On the 
other hand, when r > t* , the compaction time (t comp ), 
i.e. the time needed to reach the saturation density, is 
shorter than t, and vibrations act on the system till the 
end of the process (see Fig. 10). 




FIG. 9. Parameters to (a) and To (b) of fit Eq. (|5|) in the 
frictionless case. The squares are the numerical data, while 
the solid line (r < r*) and dotted line (r > r*) in (a) are 
the functions given in Eq. (^) and those in (b) are the func- 
tions given by Eq. (^). The horizontal shaking amplitude is 
P°h = 0.1. 

We have observed a change of behavior in correspon- 
dence of the critical value of the shaking time r*, the 
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FIG. 10 
the frictionless case 
p° = 0.1. 
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Compaction time t comp vs. shaking time r in 
The horizontal shaking amplitude is 
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We believe this is responsible for the power law tail for 
large values of r, since the combined action of shaking 
and gravity produces a slowing down of the dynamics, en- 
abling at the same time to obtain denser packings. How- 
ever, we are not sure whether this picture corresponds 
to a real dynamical phase transition in the system, and 
moreover the identification of an order parameter would 
be quite problematic. Therefore a deeper investigation 
on this point is needed and we hope to stimulate some 
future experimental and theoretical works. 

Now we briefly mention the main changes occurring 
when a sequential tapping is applied to the system. Al- 
though the single tap dynamics can give insight on the 
dynamical evolution and the relaxation of the system, 
the sequential tapping procedure can be closely related 
to experimental situations. In Fig. ^(a) we compare the 
density evolution obtained by a single tap of duration 
t = 10000 with that obtained by a sequential tapping 
applied for the same amount of time r. In particular we 
have generated 5000 taps, each of which was applied for 
a time t/5000. 



was also found in the experiments ||: 

p{t)=P -- C l + dTn{l + t/TY (10) 

In this way we can connect directly our model to the ex- 
perimental data and we can also examine the effect of dif- 
ferent dynamical shaking procedures on the compaction 
process. In this respect, our results are in agreement with 
those of Nicodemi et al. jl9| , who find a similar change of 
behavior when considering single and multiple tapping. 
This should clarify to a certain extent the mechanisms of 
density relaxation in different dynamical situations. 



IV. COMPACTION DYNAMICS WITH 
FRICTION 




FIG. 11. Comparison between single and multiple tapping. 
In (a) the behavior of density as a function of time is plot- 
ted: the solid line is obtained by a single tap with r = 10000, 
while the dashed line is the result of 5000 taps, each of which 
applied for a time r/5000. In (b) the density is shown as a 
function of the number of taps (ntap)' the solid line is the 
result of the simulation and the dashed line is a fit according 
to Eq. [k] for n tap > 40, with parameters c ~ 47, d ~ 1.2 • 10 5 
and T ~ 3 TO 4 . The horizontal shaking amplitude is p\ — 0.1. 

From the comparison we deduce that the multiple tap- 
ping dynamics is slower, suggesting that the fitting form 
Eq. [5] cannot describe properly the density evolution in 
this case. Besides, the saturation density is also lower. 
In Fig. |ll|(b) we characterize more quantitatively the dy- 
namical behavior by plotting the density as a function of 
the number of taps. It turns out that in this case the 
intermediate-long time evolution of the density can be 
described by the claimed inverse logarithmic law, which 



When a frictional force is introduced, as described in 
Sec. |, some qualitative and quantitative changes are 
found. In Fig. R2 the time evolution of the density for a 
fixed value of the friction coefficient fi and several choices 
of r is displayed. 




10 10 10 10 

t (MCS) 

FIG. 12. Temporal behavior of density for friction coeffi- 
cient = 0.4 and different values of the shaking time r (from 
left to right t = 30,50,80,150,400,2000,4000,10000). The 
horizontal shaking amplitude is pi — 0.1. 

A comparison with Fig. |^ shows that for a same value 
of r a much smaller asymptotic density is reached, and 
that a slowing down of the dynamical behavior occurs. 
As a consequence, larger values of r are needed to obtain 
the same final densities as for jj, = 0. In Fig. ^3| a com- 
parison of p(t) between /i = and \x = 0.4 for r = 80 is 
illustrated. 
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FIG. 13. Density relaxation for p = (top solid line) and 
p = 0.4 (bottom solid line). The dashed lines are stretched 
exponential fits using Eq. (H) for p(t) > 0.6, with parameters 
C ~ 1.1, t ^ 275, r ~ 400, P ~ 9 for p = 0.4 (the param- 
eters for p = can be found in the caption of Fig. pf). Wc 
note that the major differences between the two curves are 
due to the different values of j3 and poo (poo{p = 0) — 0.782 
and Poc(p ~ 0.4) ~ 0.685). The horizontal shaking amplitude 
is p° =0.1 and the shaking time is r = 80. 

It seems that the same functional form Eq. (|J) adopted 
for p = can be used to describe the dynamical evolu- 
tion of the system. In spite of that, for p ^ some slight 
discrepancies in the tails can be detected for small values 
of r, not visible at the level of detail of the figure. This 
means that the reduced probability for horizontal moves 
affects the compaction mainly in the late stages of the 
process. 

We have also examined the p dependence of p(t) . The 
curve for p, = is far above the others, and further- 
more the difference between the curves for low values of 
p is more marked. Such differences show up only after 
an initial transient in which the various curves are prac- 
tically superimposed one on the others (the curves do 
not start diverging significantly until the density reaches 
about 0.6). This is in compliance with what we have just 
said about the effectiveness of friction in the asymptotic 
dynamics. 

The behavior of the asymptotic density is depicted in 
Fig. [lJI and reveals more interesting aspects. In Fig. |l4|(a) 
we have plotted Poo(t) for p = and p = 0.4. We were 
able to identify a power law regime for small r also for 
p = 0.4: 



Poo(r) = (r/rr) 



(11) 



with ~ 0.062 (not very different from 6 in Eq. (0)) 

and t/ ~ 6.9 • 10 4 . Thus, saturation density scales in 
the same way as for p — 0, with nearly the same expo- 
nent but with a very different relaxation time (more than 
one order of magnitude larger), i.e. the process is much 



more sluggish. We did not manage to find the functional 
behavior for larger r and it is not very clear whether it 
is possible to define a crossover shaking time, at least for 
the values of r we have considered. For r > 1000 
increases more steeply than for p — and it is likely that 
one could find another regime for larger t. However com- 
putational limitations prevented us to explore in detail 
the region for r > 10000 (since simulations become more 
and more time consuming as r increases), and we will 
address this problem in the future. 
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FIG. 14. Saturation density p^ as a function of (a) r for 
p — 0.4 and (b) of p for r = 150. The circles in (a) and 
the squares in (a) and (b) are the numerical data for p — 
and for p — 0.4 respectively. The dashed lines in (a) and (b) 
are fits respectively according to Eq. ( |ll| ) and Eq. (|12|). The 
horizontal shaking amplitude is p\ — 0.1. 



Finally in Fig. |T|(b) we have analyzed the behavior of 
Poo vs. p. Interestingly poo decreases linearly with p for 
p > 0.4: 



Poo(m) = A - Bp, 



(12) 



where A ~ 0.71 and B ~ 0.073. Instead the decrease is 
more drastic for lower values of p, as noted above. This 
fact might indicate a change of behavior in the dynamics 
of the system as a function of the friction coefficient. For 
this purpose we show in Fig. |l5| the stretched exponent 
P of Eq. (|) vs. p. 
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FIG. 15. Stretched exponent j3 of Eq. (H) vs. friction coef- 
ficient fj, for pi =0.1 and r = 150. 

It rapidly decreases for fi < 0.4, but it is almost con- 
stant for larger [i. In other terms, the independence of \i 
of the dynamical evolution reflects itself just in a linear 
scaling of the final density with fj,. We hope to elucidate 
this point more thoroughly in a future work. 



totic density as a function of the shaking time with a 
slowing down above a critical shaking time value, in- 
dicating a possible dynamical transition process. More 
detailed studies of this point are planned for the future 
in order to reach a better understanding of the physical 
mechanism underlying this peculiar behavior. 

Similar to what is found in experiments, horizontal 
shaking favors locally ordered configurations and leads 
to higher compaction (2^]2^] . We believe it would be im- 
portant to address further this issue by improving our 
model, allowing rotations of the dimers during the dy- 
namics. 

The results obtained by introducing the frictional ef- 
fect in our model indicate the important role of friction 
on the compaction dynamics. In particular we find a 
slowing down of the dynamical behavior which is more 
evident in the asymptotic regime. The relaxation of the 
density is still described by a stretched exponential in the 
presence of friction, with a stretched exponent that be- 
comes constant at large values of the friction coefficient. 
It would be worthwhile to explore this issue in future 
works with an improved modeling of the frictional force. 



V. CONCLUSIONS AND PERSPECTIVES 

We have proposed and discussed a simplified model for 
dimer compaction, with the intent to take into account 
relevant dynamical features, such as horizontal shaking 
and friction forces between particles, analyzing in detail 
their effects on the density time evolution. The dynami- 
cal behavior of our model is very complex with interesting 
and novel features. In the absence of friction, for a sin- 
gle tap, the compaction dynamics cannot be interpreted 
as an inverse logarithmic relaxation but in terms of a 
stretched exponential law, which is a peculiar character- 
istic of a continuous range of time constants |3|. Thus, 
the fact that we do not find a logarithmic law but a 
stretched exponential relaxation could be a peculiar fea- 
ture related to the single tap dynamics, since Nicodemi 
et al. jl9| observe the same kind of behavior in a dif- 
ferent model with a similar dynamics. Furthermore, we 
are not aware of experiments in which a single tap is ap- 
plied as external perturbation, since they all refer to a 
tapping sequence. Actually, it turns out that the loga- 
rithmic behavior found in experiments is recovered also 
in our model when a sequential tapping is applied. We 
therefore might infer that different shaking procedures 
give rise to intrinsically different compaction behaviors, 
driven by different dynamical mechanisms. However, it 
still needs to be clarified whether the logarithmic law is 
just a good fit to the experimental data or it has a deeper 
meaning (see for example the discussion in Ref. p"^]). 

We also find a crossover in the behavior of the asymp- 



[1] Granular Matter, A. Metha ed. (Springer- Verlag 1994). 

[2] Proceedings of the NATO Advanced Study Institute on 
Physics of Dry Granular Media, Eds. H.J. Herrmann 
et al., Kluwer Academic Publishers, The Netherlands 
(1998). 

[3] H.M. Jaeger, S.R. Nagel and R.P. Behringer, Rev. Mod. 

Phys. 68, 1259 (1996). 
[4] H.M. Jaeger and S.R. Nagel, Science 255, 1523 (1992). 
[5] J.B. Knight, C.G. Fandrich, C.N. Lau, H.M. Jaeger and 

S.R. Nagel, Phys. Rev. E 51, 3957 (1995). 
[6] ER. Nowak, J.B. Knight, E. Ben-Nairn, H.M. Jaeger and 

S.R. Nagel, Phys. Rev. E 57, 1971 (1998). 
[7] ER. Nowak, J.B. Knight, M.L. Povinelli, H.M. Jaeger 

and S.R. Nagel, Powder Technol. 94, 79 (1997). 
[8] E. Ben-Naim, J.B. Knight, ER. Nowak, H.M. Jaeger and 

S.R. Nagel, Physica D 123, 380 (1998). 
[9] J. Talbot, G. Tarjus and P. Viot, Phys. Rev. E 61, 5429 

(2000) . 

[10] S.J. Linz, Phys. Rev. E 54, 2925 (1996). 
[11] S.F. Edwards and D.V. Grinev, Phys. Rev. E 58, 4758 
(1998). 

[12] J.J. Brey, A. Prados and B. Sanchez-Rey, Phys. Rev. E 

60, 5685 (1999). 
[13] K.L. Gavrilov, Phys. Rev. E 58, 2107 (1998). 
[14] P. Philippe and D. Bideau, Phys. Rev. E 63, 051304 

(2001) . 

[15] D.A. Head and G.J. Rodgers, J. Phys. A 31, 107 (1998). 

[16] D.A. Head, Phys. Rev. E 62, 2439 (2000). 

[17] S. Luding, M. Nicolas and O. Pouliquen, in Compaction 
of soils, Granulates and Powders, D. Kolymbas, W. Fellin 
and A. A. Balkema (eds.), p. 241, Rotterdam (2000). 



9 



[18] J.-M. Hertzsch, Eur. Phys. J. B 18, 459 (2000). 

[19] M. Nicodemi, A. Coniglio and H.J. Herrmann, Phys. Rev. 

E 55, 3962 (1997). 
[20] M. Nicodemi, Phys. Rev. Lett. 82, 3734 (1999). 
[21] M. Nicodemi, Physica A 285, 267 (2000). 
[22] E. Caglioti, V. Loreto, H.J. Herrmann and M. Nicodemi, 

Phys. Rev. Lett. 79, 1575 (1997). 
[23] A. Barrat and V. Loreto, J. Phys. A 33, 4401 (2000). 
[24] A. Barrat, J. Kurchan, V. Loreto and M. Sellitto, Phys. 

Rev. Lett. 85, 5034 (2000). 
[25] M. Sellitto and J.J. Arenzon, Phys. Rev. E 62, 7793 

(2000). 

[26] W. Kob and H. C. Andersen, Phys. Rev. E 48, 4364 
(1993). 

[27] Y. Levin , J.J. Arenzon and M. Sellitto, Europhys. Lett. 



55, 767 (2001). 

[28] O. Pouliquen, M. Nicolas and P.D. Weidman, Phys. Rev. 

Lett. 79, 3640 (1997). 
[29] M. Nicolas, P. Duru and O. Pouliquen, Eur. Phys. J. E 

3, 309 (2000). 

[30] F.X. Villarruel, B.E. Lauderdale, D.M. Mueth and H.M. 

Jaeger, Phys. Rev. E 61, 6914 (2000). 
[31] C. Fusco, P. Gallo, A. Petri and M. Rovere, Phys. Rev. 

E 65, 026127 (2002). 
[32] P.W. Kasteleyn, Physica 27, 1209 (1961). 
[33] M.E. Fisher, Phys. Rev. 124, 1664 (1961). 
[34] C. Fusco, P. Gallo, A. Petri and M. Rovere, J. Chem. 

Phys. 114, 7563 (2001). 
[35] J. Duran, T. Mazozi, E. Clement and J. Rajchenbach, 

Phys. Rev. E 50, 3092 (1994). 



10 



